clear
rng default
addpath /Users/cristinagualdani/mosek/9.2/toolbox/r2015aom
%% Load useful matrices
load PXY_cond
load V_CS


%% Construct grid of parameters to check for each realisation of X
u0grid_temp=0; %[Ychosen=0, X=1] [Ychosen=0, X=2], location normalisation 
u1grid_temp_all=[V_CS(1,1) V_CS(1,2) V_CS(1,3)]; %[Ychosen=1, X=1] [Ychosen=1, X=2] [Ychosen=1, X=3], scale normalisation
gran=0.1;  
u2grid_temp=-20:gran:20;  %[Ychosen=2, X=1] [Ychosen=2, X=2] [Ychosen=2, X=3]
u3grid_temp=-20:gran:20;  %[Ychosen=3, X=1] [Ychosen=3, X=2] [Ychosen=3, X=3]
% try different values of gran and different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for higher values of gran and larger sizes of hypercube)

%% Construct identified set for each realisation of X
ntypes_w=3;
IdSetW=cell(ntypes_w,1);


for x=1:ntypes_w %we can proceed separately across types as long as we do not impose independence from X
    u1grid_temp=u1grid_temp_all(x);
    [ca, cb, cc,cd] = ndgrid(u0grid_temp, u1grid_temp, u2grid_temp, u3grid_temp);
    u0grid=ca(:);
    u1grid=cb(:);
    u2grid=cc(:);
    u3grid=cd(:);
    sg=size(u0grid,1);
    % Components invariant across parameter values
    [sgreduced, ...
    id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
    P_col,...
    Tcoord_zeros,...
    Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
    cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros,...
    equalorder, n_vertices,infty_pos,...
    cr_symmetry, fill_symmetry, RHS_symmetry, S1, S2, S3]=useful_anyparam(u0grid, u1grid, u2grid, u3grid, sg);
    % Loop over parameter values, order of columns: [u0x u1x u2x u3x]
    IdSetW{x}=OneSide(u0grid, u1grid, u2grid, u3grid, ...
                       sgreduced, PXY_cond(x,:), ...
                       id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
                       P_col,...
                       Tcoord_zeros,...
                       Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
                       cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros, equalorder, n_vertices,infty_pos,...
                       cr_symmetry, fill_symmetry, RHS_symmetry, S1, S2, S3);  
end

save('IdSetW.mat', 'IdSetW')


